/*!  \file  Matrix.h
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      04/21/2013 05:17:33 PM
 *   \copyright GNU Public License.
 */
#pragma once


// An m-by-n Matrix has m rows and n columns
// For example the following 2-by-3 matrix
//  1 2 3
//  4 5 6
//  is stored as {1 2 3 4 5 6} in RowMajor format
//  and as       {1 4 2 5 3 6} in ColMajor format
//  the opperation of changing the data layout is
//  called in-place transpose of a Matrix

// This class is designed to have full capability of hangling all kind of
// situations where matrix/matrix operations are needed. However, I will only
// implement the once that I need for my master thesis. If someone else uses 
// this code, make sure that the operation you are using is implementeted.
typedef unsigned int uint;
typedef const unsigned int cuint; // short cut
#include "macros.h"
class Matrix {
public:
	enum Type 		  { General = 0, Orthogonal = 1, Diagonal = 2, BiDiagonal = 3, TriDiagonal = 4, GivensRotations = 5, ElementaryReflectors = 6};
	// Row-major and column-major order describe methods for  storing arrays in
	// linear memory. C uses RowMajor
	enum Order		  { ColMajor, RowMajor};
	// Native := no internal structure
	// ElementaryReflectors := Elementary Householder matrices 
	enum Operation   { Transpose, TransposeData, Copy, CopyWithData, Unknown };
// CONSTRUCTORS	
	Matrix(const int m_, const int n_, const Type type_ = General, 
				const double * const data_ = 0);
	Matrix(const Matrix & m);
	Matrix(const Matrix & m, const Operation  op);
	virtual ~Matrix(){ delocateStorage(); };
	bool allocateStorage();
	void delocateStorage();

	void transpose    ();						 // transposes the matrix O(1)
	bool changeOrder  ();						 // changes the actual memory layout
	void copy   		(const Matrix & mat); // copy the entire matrix. 
	void copyWithData (const Matrix & mat); // new storage and copy data 
	void copyOnlyData (const Matrix & mat); // keep the current storage
	Matrix T() const { return Matrix (*this, Transpose); }
	void initialize(const double * const source);
	void fill(const double val);

	Matrix & operator  = (const Matrix & m);
	Matrix & operator *= (const Matrix & m);
	Matrix & operator += (const Matrix & m) {add(m, 1.0); return *this; }
	Matrix & operator -= (const Matrix & m) {add(m,-1.0); return *this; }

	void add(const Matrix & mat, const double alpha = 1.0);
	bool makeGeneral(){ type = isDense() ? General : type; return isGeneral();}
	bool makeOrthogonal(){ type = isDense() ? Orthogonal : type; return isOrthogonal();}
	void makeRowMajor();
	void makeColMajor();
	void makeUpper()	  { makeRowMajor();   };	
	void makeLower()	  { makeColMajor();   };	
	void makeIdentity();
	void makeQ()		  { assert(isReflector()); Q = true; };
	void makeP()		  { assert(isReflector()); Q = false; };
	void setLD()        { ld = isColMajor() ? M : N;};
	void setRefl(const Matrix & m, const bool Q1) {
		assert(isReflector());
		relflectorStorage = &m; 
		Q = Q1;
	};

	void generateOrthogonal(); // reflector to orthogonal
 
	double 		 * getModData()    { 
		if (count() > 1 && isFullMatrix() ) {
			 *this = Matrix(*this, CopyWithData); }
		return data; }
	double       * getModOffDiagonal() { 
						return getModData() + std::min(N,M);	}

// Matrix manipulations
	friend std::ostream& operator<<(std::ostream& os, const Matrix & m);
	const Matrix operator*(const Matrix & m) const;
	const Matrix operator+(const Matrix & m) const {
		return (Matrix(*this,CopyWithData) += m);};
	const Matrix operator-(const Matrix & m) const {
		return (Matrix(*this,CopyWithData) -=m);};

	void applyGivensRotationFromLeft (const unsigned int i, 
		const unsigned int j, const double c, const double s);

	void applyGivensRotationFromRight(const unsigned int i, 
		const unsigned int j, const double c, const double s);

	void killEntry(const unsigned int i, const unsigned int j, 
				 const unsigned int q, bool fromLeft);

	void flipCols();// columns flipped in the left-right direction
	void flipRows();// rows    flipped in the up-down direction
	void changeColSign(const unsigned int i);

	// Sub Matrix setting
	void setView(const int i_, const int j_, const int m_, const int n_);
	void resetView();
	

	void deflateLeft (Matrix & L); // deflates the matrix from left
	void deflateRight(Matrix & R); // deflates the matrix from right
	void deflate     (Matrix & L, Matrix & R); // deflates from left & right

	void setRotations(const int maxNumberOfRotations);

	void applyReflector(const Matrix & Refl, const bool fromLeft);
	void applyReflectorFromLeft (const Matrix & Refl) {
		applyReflector(Refl, true);
	}
	void applyReflectorFromRight(const Matrix & Refl) {
		applyReflector(Refl, false);
	}

	bool sameDim    (const Matrix& m) const { return sameRows(m)&&sameCols(m);}
	bool sameRows	 (const Matrix& m) const { return rows()  == m.rows() ; }
	bool sameCols	 (const Matrix& m) const { return cols()  == m.cols() ; }
	bool sameOrder  (const Matrix& m) const {	return order   == m.order  ; }
	bool sameType   (const Matrix& m) const {	return type    == m.type   ; }
	bool sameDense  (const Matrix& m) const {	return isDense() == m.isDense(); }
	bool same       (const Matrix& m) const { return sameType(m) && sameDim(m) && sameOrder(m);	}

	bool isDense()      const { return isGeneral() || isOrthogonal() || 
												  isReflector(); };
	bool isSparse() 	  const { return !isDense(); };
	bool isGeneral()    const { return type  == General    ; };
	bool isOrthogonal() const { return type  == Orthogonal ; };
	bool isReflector () const { return type  == ElementaryReflectors; }; 
	bool isDiagonal()   const { return type  == Diagonal   ; };
	bool isBiDiagonal() const { return type  == BiDiagonal ; };
	bool isTriDiagonal()const { return type  == TriDiagonal; };
	bool isRowMajor()   const { return order == RowMajor; }
	bool isColMajor()   const { return order == ColMajor; }
	bool isUpper()		  const { return isRowMajor(); };	
	bool isLower()		  const { return isColMajor(); };	
	bool isSquare()	  const { return rows() == cols(); };	
	bool isSubMatrix()  const { return !isFullMatrix(); };
	bool isFullMatrix() const { return n == N && m == M;};
	bool isQ()		     const { assert(isReflector()); return  Q; };
	bool isP()		     const { assert(isReflector()); return !Q; };
	Type getType()      const { return type; };
	Order getOrder()    const { return order; };
	

	//  The entry of a matrix A is written using two indices, A_i,j
	//  where the first subscript is the row number 
	//  and the second is the column number.
	uint offsetColMajor(cuint row, cuint col)	const {return col * ld + row;}
	uint offsetRowMajor(cuint row, cuint col) const {return row * ld + col;}
	uint offset        (cuint row, cuint col) const {return isRowMajor() ? 
			offsetRowMajor(row, col) : offsetColMajor(row, col);}
	void getRowColFromOffset(const int offset, int & row, int & col) const {
		div_t q = div(offset, ld);
		row = isRowMajor() ? q.quot : q.rem;
		col = isColMajor() ? q.quot : q.rem;
	}
	const double & operator() (const int row, const int col) const {
		return getData()[offset(row,col)];	}; 

	int rows() const { return m; };
	int cols() const { return n; };
	// leading dimention := distance in memory between elements of two
	// consecutive columns which have the same row index.
	int lda () const { assert (ld == (isRowMajor()? N : M));return ld;}
	int gap () const { return ld - (isRowMajor()? cols() : rows()); }
	int minDim() const { return n < m ? n : m; }
	uint getStorageSize() const;
	int count () const { return counter ? *counter : -1; };	

	// Entrywise norms
	double 		nrm2() const;// returns Frobenius norm 
	double  OrthNorm() const;// max_{i!=j} |v_i . v_j|

	const double * getData() const { return data; }
	const double * getOffDiagonal   () const { 
						return data         + std::min(N,M); }


	void solve(const double * RS b, const int ldb, double * RS x, 
		const int ldx,const double alpha = 1.0) const; // solves A x = b

	 
	const Matrix * getRefl() const { return relflectorStorage; };
// class members
// for the future make all member private
private:
	void swapCols(const unsigned int i, const unsigned int j);
	void swapRows(const unsigned int i, const unsigned int j);
	void DiagMultCols(const double * diag); // col[i] *= diag[i]
	void DiagMultRows(const double * diag); // row[i] *= diag[i]
	void rotate(const unsigned int nSteps, const unsigned int step, 
			const double c, const double s, double * RS i, double * RS j);
	
// data members
	double *    data;
	int	 *    counter;
	int         m; 			// number of rows in current view
	int         n; 			// number of cols in current view
	int 			ld; 			// leading dimetion of matrix
	int 			shift;      // >= 0 for subMatrix
	int			M;				// number of rows in full matrix
	int			N;				// number of cols in full matrix
	Order			order;
	Type			type;
	const Matrix * relflectorStorage; // != NULL, if the matrix is reflector 
	int 		    * rotationIndex; // used with Givens rotatation in packed form
	bool 			Q;				// if the matrix is reflecto is it Q or P
};
const Matrix::Type     General          = Matrix::General;
const Matrix::Type    Diagonal          = Matrix::Diagonal;
const Matrix::Type  Orthogonal          = Matrix::Orthogonal;
const Matrix::Type  BiDiagonal          = Matrix::BiDiagonal;
const Matrix::Type TriDiagonal          = Matrix::TriDiagonal;
const Matrix::Type GivensRotations      = Matrix::GivensRotations;
const Matrix::Type ElementaryReflectors = Matrix::ElementaryReflectors;

const Matrix::Order RowMajor   			 = Matrix::RowMajor;
const Matrix::Order ColMajor   			 = Matrix::ColMajor;

const Matrix::Operation Transpose     = Matrix::Transpose    ;
const Matrix::Operation TransposeData = Matrix::TransposeData;
const Matrix::Operation Copy          = Matrix::Copy         ;
const Matrix::Operation CopyWithData  = Matrix::CopyWithData ;

#include <math.h>
inline void rot(double & a, double & b, const double & c, const double & s) {
	double tmp1 = a * c  + b * s;
	double tmp2 = a * -s + b * c;
	a = tmp1; b = tmp2;
}
inline void getCS(const double & a, const double & b, double & c, double & s) {
	if (fabs(b) > fabs(a)) {
		double t = a/b;
 		double u = copysign(sqrt(1+t*t),b);
		s = 1/u;
  		c = s*t;
	} else {
  		double t =  b/a;
  		double u = copysign(sqrt(1+t*t),a);
  		c = 1/u;
  		s = c*t;	
	}
}
inline void getCS2(const double & a, const double & b, double & c, double & s) {
	bool cond = fabs(b) > fabs(a);
	double in[2]; 
	in[cond] = a; in[!cond] = b;
	double t = in[1] / in[0];
 	double u = copysign(sqrt(1+t*t),in[!cond]);
	in[0] = 1.0/u;
	in[1] = in[0] * t;
	c = in[!cond];
	s = in[ cond];
}
